suppressPackageStartupMessages({
library(Seurat, quietly = T)
library(monocle3, quietly = T)
library(SeuratWrappers, quietly = T)
library(WGCNA, quietly = T)
library(hdWGCNA, quietly = T)
library(clusterProfiler)
library(dplyr, quietly = T)
library(tidyr, quietly = T)
library(magrittr, quietly = T)
library(tibble, quietly = T)
library(reshape2, quietly = T)
library(stringr)
library(plyr, quietly = T)
library(textshape, quietly = T)
library(ggplot2, quietly = T)
library(cowplot, quietly = T)
library(RColorBrewer, quietly = T)
library(patchwork)
})
server = 'mando'
if (server == 'jabba'){
data_path = '/data3/hratch/norcross_abc/'
}else if (server == 'mando'){
data_path = '/data/hratch/norcross_abc/'
}
n.cores<-20
seed<-888
set.seed(12345)
# optionally enable multithreading
WGCNA::enableWGCNAThreads(nThreads = n.cores)
# Load and format T cell data to get just the CD8s
abc.tcells<-readRDS(paste0(data_path, 'processed/abc_tcells.RDS'))
Idents(abc.tcells)<-'Cell.Type.Level2'
abc.tcells
An object of class Seurat 20723 features across 14214 samples within 2 assays Active assay: integrated (2000 features, 2000 variable features) 1 other assay present: RNA 2 dimensional reductions calculated: pca, umap
Subset to just the CD8+ T cells
abc.cd8s.all<-abc.tcells[, unname(sapply(as.character(abc.tcells$Cell.Type.Level2), function(x) startsWith(x, 'CD8+')))]
We want to place the root node at the CD8 naive T-cells Looking at the UMAP, we can see there is a main group of CD8+ naive T-cells (clusters 30, 41, and 43) and then an additional cluster (cluster 44) that is transcriptionally distinct in UMAP space. Accordingly, we would like to set the root node based on the main cluster:
h_ = 5
w_ = 12
options(repr.plot.height=h_, repr.plot.width=w_)
g1A<-DimPlot(abc.cd8s.all)
g1B<-DimPlot(abc.cd8s.all, group.by = 'Seurat.Clusters.Level2')
g1<-cowplot::plot_grid(g1A, g1B, ncol = 2)
g1
h_ = 10
w_ = 20
options(repr.plot.height=h_, repr.plot.width=w_)
g2A<-DimPlot(abc.cd8s.all, split.by = 'orig.ident')
g2B<-DimPlot(abc.cd8s.all, group.by = 'Seurat.Clusters.Level2', split.by = 'orig.ident')
g2<-cowplot::plot_grid(g2A, g2B, ncol = 1)
g2
md<-abc.cd8s.all@meta.data
md[['Root.Node']]<-'no'
md[(md$Cell.Type.Level2 == 'CD8+ TN') & (md$Seurat.Clusters.Level2 != '44'), 'Root.Node']<-'yes'
abc.cd8s.all@meta.data<-md
Do the trajectory analysis:
# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, root.node="yes"){
cell_ids <- which(colData(cds)[, "Root.Node"] == root.node)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
Run trajectory analysis on all samples simultaneously rather than each one separately:
abc.cd8s.all.cds<-SeuratWrappers::as.cell_data_set(abc.cd8s.all)
abc.cd8s.all.cds<-monocle3::cluster_cells(cds = abc.cd8s.all.cds)
abc.cd8s.all.cds<-monocle3::learn_graph(cds = abc.cd8s.all.cds)#, use_partition = F)
abc.cd8s.all.cds <- order_cells(abc.cd8s.all.cds, root_pr_nodes=get_earliest_principal_node(abc.cd8s.all.cds))
saveRDS(abc.cd8s.all.cds, paste0(data_path, 'processed/', 'cd8s_monocole_obj.rds'))
|======================================================================| 100%
saveRDS(abc.cd8s.all.cds, paste0(data_path, 'processed/', 'cd8s_monocole_obj.rds'))
white circles are roots, gray circles are leaves, blue circles are branch points:
h_ = 6
w_ = 12
options(repr.plot.height=h_, repr.plot.width=w_)
g1A<-DimPlot(abc.cd8s.all) + ggtitle('All Conditions')
g1B<-plot_cells(cds = abc.cd8s.all.cds, color_cells_by = "pseudotime", show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE)
g1<-g1A+g1B
for (ext in c('.svg', '.png', '.pdf')){
fn<-paste0(data_path, 'figures/', 'cd8_trajectory_all', ext)
ggsave(fn, g1, height = h_, width = w_)}
g1
traj.plots<-list()
conditions <- unique(abc.cd8s.all.cds$orig.ident)
for (i in seq_along(conditions)){
g1<-DimPlot(abc.cd8s.so[[i]]) + ggtitle(paste0('Condition: ', conditions[[i]]))
# if (i != length(abc.cd8s.all.cds)){
# g1<-g1 + theme(legend.position="none")
# }
g2<-plot_cells(cds = abc.cd8s.all.cds[, abc.cd8s.all.cds$orig.ident == conditions[[i]]],
color_cells_by = "pseudotime", show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE)
g<-cowplot::plot_grid(g1, g2, ncol = 1)#g1 + g2
traj.plots[[conditions[[i]]]]<-g
}
h_ = 12
w_ = 35
options(repr.plot.height=h_, repr.plot.width=w_)
g2<-cowplot::plot_grid(plotlist=traj.plots, ncol = length(traj.plots))
for (ext in c('.svg', '.png', '.pdf')){
fn<-paste0(data_path, 'figures/', 'cd8_trajectory_all_bycondition', ext)
ggsave(fn, g2, height = h_, width = w_)}
g2
from: https://cole-trapnell-lab.github.io/monocle3/docs/differential/#branches
For branches points, we take cells a little before (psuedotime-wise) the branch point, and cells from each path a little after the branch point.
For paths between nodes, we just take the cells along the trajectory immediately after node 1 (assuming node 1 is a branch) to immediately before node 2 if node 2 is a branch, or cells radially sorrounding node 2 if node 2 is a terminal state.
These regions are selected visually on the UMAP
abc.cd8s.all.cds <- readRDS(paste0(data_path, 'processed/', 'cd8s_monocole_obj.rds'))
abc.cd8s.all.cds <- estimate_size_factors(abc.cd8s.all.cds) # size factors needed for branching analysis
get_subsets<-function(cds, branch_umap, branch_umap_exclude = NULL){
umap.coords<-reducedDims(abc.cd8s.all.cds)[['UMAP']]
if (!(is.null(branch_umap_exclude))){
cell.ids.exclude<-c()
for (bue in branch_umap_exclude){
cond1 <- (umap.coords[, 1] >= bue$UMAP_1[1])
cond2 <- (umap.coords[, 1] <= bue$UMAP_1[2])
cond3 <- (umap.coords[, 2] >= bue$UMAP_2[1])
cond4 <- (umap.coords[, 2] <= bue$UMAP_2[2])
cell.ids.exclude<-c(cell.ids.exclude,
rownames(umap.coords[cond1 & cond2 & cond3 & cond4,]))
}
}
cond1 <- (umap.coords[, 1] >= branch_umap$UMAP_1[1])
cond2 <- (umap.coords[, 1] <= branch_umap$UMAP_1[2])
cond3 <- (umap.coords[, 2] >= branch_umap$UMAP_2[1])
cond4 <- (umap.coords[, 2] <= branch_umap$UMAP_2[2])
if (!(is.null(branch_umap_exclude))){
cond5 <- (!(rownames(umap.coords) %in% cell.ids.exclude))
all.conds <- cond1 & cond2 & cond3 & cond4 & cond5
}else{
all.conds <- cond1 & cond2 & cond3 & cond4
}
cell.ids<-rownames(umap.coords[all.conds,])
return(cell.ids)
}
We will visually select cell subsets according to UMAP coordinates to analyze branches/trajectories of interest:
branch.cells.list<-list()
Branch 1:
branch.label<-'Branch1.cells'
branch_umap<-list('UMAP_1' = c(-4.25, -2), 'UMAP_2' = c(-2, 0.3))
branch_umap_exclude<-list(list('UMAP_1' = c(-4.25, -3.5), 'UMAP_2' = c(-2, -1)),
list('UMAP_1' = c(-3.5, -2), 'UMAP_2' = c(-2, -1.5)))
branch.cells<-get_subsets(abc.cd8s.all.cds, branch_umap, branch_umap_exclude)
colData(abc.cd8s.all.cds)[[branch.label]]<-'no'
colData(abc.cd8s.all.cds)[branch.cells, branch.label]<-'yes'
colData(abc.cd8s.all.cds)[[branch.label]] <- factor(colData(abc.cd8s.all.cds)[[branch.label]], levels = c('yes', 'no'))
branch.cells.list[[branch.label]]<-branch.cells
h_ = 6
w_ = 6
options(repr.plot.height=h_, repr.plot.width=w_)
g.branch<-monocle3::plot_cells(cds = abc.cd8s.all.cds, color_cells_by = branch.label, show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE) + ggtitle(branch.label)
g.branch
Branch 3
branch.label<-'Branch3.cells'
branch_umap<-list('UMAP_1' = c(-2.25, 1), 'UMAP_2' = c(-1.75, 0.5))
branch_umap_exclude<-NULL#list(list('UMAP_1' = c(-2.75, -1.75), 'UMAP_2' = c(-2.25, -0.25)))
branch.cells<-get_subsets(abc.cd8s.all.cds, branch_umap, branch_umap_exclude)
colData(abc.cd8s.all.cds)[[branch.label]]<-'no'
colData(abc.cd8s.all.cds)[branch.cells, branch.label]<-'yes'
colData(abc.cd8s.all.cds)[[branch.label]] <- factor(colData(abc.cd8s.all.cds)[[branch.label]], levels = c('yes', 'no'))
branch.cells.list[[branch.label]]<-branch.cells
h_ = 6
w_ = 6
options(repr.plot.height=h_, repr.plot.width=w_)
g.branch<-plot_cells(cds = abc.cd8s.all.cds, color_cells_by = branch.label, show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE) + ggtitle(branch.label)
g.branch
Path 1A: from branch 3 to terminal state 3
branch.label<-'Path1A.cells'
branch_umap<-list('UMAP_1' = c(-1.25, 1.5), 'UMAP_2' = c(-4, -1))
branch_umap_exclude<-NULL#list(list('UMAP_1' = c(-2.75, -1.75), 'UMAP_2' = c(-2.25, -0.25)))
branch.cells<-get_subsets(abc.cd8s.all.cds, branch_umap, branch_umap_exclude)
colData(abc.cd8s.all.cds)[[branch.label]]<-'no'
colData(abc.cd8s.all.cds)[branch.cells, branch.label]<-'yes'
colData(abc.cd8s.all.cds)[[branch.label]] <- factor(colData(abc.cd8s.all.cds)[[branch.label]], levels = c('yes', 'no'))
branch.cells.list[[branch.label]]<-branch.cells
h_ = 6
w_ = 6
options(repr.plot.height=h_, repr.plot.width=w_)
g.branch<-plot_cells(cds = abc.cd8s.all.cds, color_cells_by = branch.label, show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE) + ggtitle(branch.label)
g.branch
Path 1B: from branch 3 to branch 4
branch.label<-'Path1B.cells'
branch_umap<-list('UMAP_1' = c(-2.25, 0), 'UMAP_2' = c(-1, 2.25))
branch_umap_exclude<-list(list('UMAP_1' = c(-2.25, -1.75), 'UMAP_2' = c(-1, -0.25)),
list('UMAP_1' = c(-2.25, -1.25), 'UMAP_2' = c(1.75, 2.25)))
branch.cells<-get_subsets(abc.cd8s.all.cds, branch_umap, branch_umap_exclude)
colData(abc.cd8s.all.cds)[[branch.label]]<-'no'
colData(abc.cd8s.all.cds)[branch.cells, branch.label]<-'yes'
colData(abc.cd8s.all.cds)[[branch.label]] <- factor(colData(abc.cd8s.all.cds)[[branch.label]], levels = c('yes', 'no'))
branch.cells.list[[branch.label]]<-branch.cells
h_ = 6
w_ = 6
options(repr.plot.height=h_, repr.plot.width=w_)
g.branch<-plot_cells(cds = abc.cd8s.all.cds, color_cells_by = branch.label, show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE) + ggtitle(branch.label)
g.branch
Path 2: from branch 4 to terminal state 2
branch.label<-'Path2.cells'
branch_umap<-list('UMAP_1' = c(0.25, 10), 'UMAP_2' = c(-1.5, 5))
branch_umap_exclude<-list(list('UMAP_1' = c(0.25, 0.75), 'UMAP_2' = c(-2.25, 1.5)),
list('UMAP_1' = c(0.25, 4), 'UMAP_2' = c(-2.25, 1)))
branch.cells<-get_subsets(abc.cd8s.all.cds, branch_umap, branch_umap_exclude)
colData(abc.cd8s.all.cds)[[branch.label]]<-'no'
colData(abc.cd8s.all.cds)[branch.cells, branch.label]<-'yes'
colData(abc.cd8s.all.cds)[[branch.label]] <- factor(colData(abc.cd8s.all.cds)[[branch.label]], levels = c('yes', 'no'))
branch.cells.list[[branch.label]]<-branch.cells
h_ = 6
w_ = 6
options(repr.plot.height=h_, repr.plot.width=w_)
g.branch<-plot_cells(cds = abc.cd8s.all.cds, color_cells_by = branch.label, show_trajectory_graph = TRUE,
trajectory_graph_segment_size = 1, trajectory_graph_color = "black",
graph_label_size = 3,
label_cell_groups=FALSE, label_roots = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE) + ggtitle(branch.label)
g.branch
q.val.thresh<-0.1
metascape_input = list()
# note, CDS only uses the 2000 HVGs
metascape_input[['_BACKGROUND']]<-rownames(abc.cd8s.all.cds) # the universe
der.branches<-list()
for (branch.label in names(branch.cells.list)){
branch.cells<-branch.cells.list[[branch.label]]
# do the DE analysis (on the 2k HVGs in the integrated counts)
# doing on integrated counts is fine bc unlike transcriptional analysis script 4, there is no controlling
# for batch in this moran's I "graph_test"
abc.cd8s.branch.cds<-abc.cd8s.all.cds[, branch.cells]
der.branch <- graph_test(abc.cd8s.branch.cds, neighbor_graph="principal_graph", cores=n.cores)
der.branch <- der.branch[(der.branch$status == 'OK'),]# & (der.branch$q_value <= q.val.thresh),]
der.branch<-der.branch[with(der.branch, order(-abs(morans_I), q_value)), ]
der.branch[['Condition']]<-branch.label
der.branch[['gene.id']]<-rownames(der.branch)
der.branches[[branch.label]]<-der.branch
metascape_input[[branch.label]]<-rownames(der.branch[der.branch$q_value <= q.val.thresh,])
# de.genes<-rownames(der.branches)
# # run coexpression analysis on
# # run on all HVGs, not just the identified DE genes, otherwise too few modules are identified
# wgcna.branch<-abc.cd8s.all[,branch.cells] # subset to branch
# # remove all 0 genes
# expr<-wgcna.branch@assays$RNA@data
# genes.to.keep<-rownames(expr[rowSums(wgcna.branch[])>0, ])
# wgcna.branch<-wgcna.branch[genes.to.keep, ]
# wgcna.branch <- hdWGCNA::SetupForWGCNA(seurat_obj = wgcna.branch,
# gene_select = 'custom', gene_list = de.genes,
# # gene_select = 'variable',
# wgcna_name = branch.label)
# wgcna.branch <- hdWGCNA::SetDatExpr(seurat_obj = wgcna.branch,
# group_name = as.character(unique(abc.cd8s.branch.so$Cell.Type.Level2)),
# group.by='Cell.Type.Level2', assay = 'RNA', slot = 'data')
# wgcna.branch <- hdWGCNA::TestSoftPowers(wgcna.branch, networkType = 'signed')
# wgcna.branch <- hdWGCNA::ConstructNetwork(wgcna.branch, setDatExpr=FALSE, overwrite_tom = TRUE,
# soft_power=NULL, # automated selection
# )
# modules <- GetModules(wgcna.branch)
}
# format input to metascape
max.length <- max(sapply(metascape_input, length))
metascape_input <- lapply(metascape_input, function(v) { c(v, rep(NA, max.length-length(v)))})
metascape_input<-do.call(cbind, metascape_input)
write.csv(metascape_input, paste0(data_path, 'interim/', 'traj_analysis_metascape_input.csv'), row.names=FALSE)
|=======================================================| 100%, Elapsed 00:03 |=======================================================| 100%, Elapsed 00:03 |=======================================================| 100%, Elapsed 00:02 |=======================================================| 100%, Elapsed 00:03 |=======================================================| 100%, Elapsed 00:05
de.res<-do.call("rbind", der.branches)
de.res<-de.res[de.res$q_value <= q.val.thresh, ]
write.csv(de.res, paste0(data_path, 'processed/', 'traj_analysis_de.csv'))
No. of differentially expressed genes per condition (note, only the 2k HVGs are run in trajectory analysis):
table(de.res$Condition)
Branch1.cells Branch3.cells Path1A.cells Path1B.cells Path2.cells
94 158 145 141 68
Run metascape on the DE genes (q <= 0.1) on all conditions at once and visualize output below:
format.ms<-function(ms, val_col){
val_df<-as.data.frame(pivot_wider(data = ms, id_cols = 'GeneList', names_from = 'Description',
values_from = val_col,
values_fn = mean))
val_df<-t(column_to_rownames(x = val_df, loc = 1))
return(val_df)
}
ms.top20<-read.csv(paste0(data_path, 'interim/', 'traj_analysis_metascape_out/Enrichment_heatmap/HeatmapSelectedGO.csv'))
ms.top100<-read.csv(paste0(data_path, 'interim/', 'traj_analysis_metascape_out/Enrichment_heatmap/HeatmapSelectedGOTop100.csv'))
ms<-read.csv(paste0(data_path, 'interim/', 'traj_analysis_metascape_out/Enrichment_GO/GO_AllLists.csv'))
First, let's see if there is any duplicate Description/GeneList combinations:
ms %>%
dplyr::group_by(GeneList, Description) %>%
dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
dplyr::filter(n > 1L)
| GeneList | Description | n |
|---|---|---|
| <chr> | <chr> | <int> |
There are no duplicate terms
freq<-format.ms(ms = ms, val_col = 'X.InGO')[ms.top20$Description, ]
pvals<-format.ms(ms = ms, val_col = 'Log.q.value.')[ms.top20$Description, ]
pvals<-as.data.frame(pvals)
pvals[['Enrichment.Term']]<-rownames(pvals)
freq<-as.data.frame(freq)
freq[['Enrichment.Term']]<-rownames(freq)
pvals<-melt(pvals, id.vars = 'Enrichment.Term', value.name = 'log10p', variable.name = 'Condition')
freq<-melt(freq, id.vars = 'Enrichment.Term', value.name = 'Frequency', variable.name = 'Condition')
pvals[['log10p']] <- -as.numeric(pvals[['log10p']])
viz.df <- cbind(pvals, Frequency = as.numeric(freq$Frequency))
h_ = 10
w_ = 25
options(repr.plot.height=h_, repr.plot.width=w_)
green_hex = brewer.pal(n = 11, name ='RdYlGn')[[11]]
yellow_hex = brewer.pal(n = 11, name = 'RdYlGn')[[6]]
red_hex = brewer.pal(n = 11, name = 'RdYlGn')[[1]]
lower_q = min(viz.df$log10p[!(is.na(viz.df$log10p))])
upper_q = max(viz.df$log10p[!(is.na(viz.df$log10p))])
middle_q = mean(c(lower_q, upper_q)) # median(viz.df$log10p[!(is.na(viz.df$log10p))]) #
# middle_q = 7.5
g<-ggplot(data = viz.df, aes(x = Condition, y = Enrichment.Term, color = log10p, size = Frequency)) +
geom_point() +
# scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = '-log10(q-val)',
# limits = c(lower_q, upper_q), midpoint = middle_q) +
scale_colour_gradientn(colours = rev(c("darkred", "orange", "yellow")), name = '-log10(q-val)')+
scale_size_continuous(range = c(1,15)) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
text = element_text(size = 20))
for (ext in c('.svg', '.png', '.pdf')){
fn<-paste0(data_path, 'figures/', 'traj_analysis_metascape_enrichment_dotplot', ext)
ggsave(fn, g, height = h_, width = w_)}
g
Warning message: “Removed 59 rows containing missing values (`geom_point()`).” Warning message: “Removed 59 rows containing missing values (`geom_point()`).” Warning message: “Removed 59 rows containing missing values (`geom_point()`).” Warning message: “Removed 59 rows containing missing values (`geom_point()`).”
freq<-format.ms(ms = ms, val_col = 'X.InGO')[ms.top100$Description, ]
pvals<-format.ms(ms = ms, val_col = 'Log.q.value.')[ms.top100$Description, ]
pvals<-as.data.frame(pvals)
pvals[['Enrichment.Term']]<-rownames(pvals)
freq<-as.data.frame(freq)
freq[['Enrichment.Term']]<-rownames(freq)
pvals<-melt(pvals, id.vars = 'Enrichment.Term', value.name = 'log10p', variable.name = 'Condition')
freq<-melt(freq, id.vars = 'Enrichment.Term', value.name = 'Frequency', variable.name = 'Condition')
pvals[['log10p']] <- -as.numeric(pvals[['log10p']])
viz.df <- cbind(pvals, Frequency = as.numeric(freq$Frequency))
h_ = 20
w_ = 25
options(repr.plot.height=h_, repr.plot.width=w_)
green_hex = brewer.pal(n = 11, name ='RdYlGn')[[11]]
yellow_hex = brewer.pal(n = 11, name = 'RdYlGn')[[6]]
red_hex = brewer.pal(n = 11, name = 'RdYlGn')[[1]]
lower_q = min(viz.df$log10p[!(is.na(viz.df$log10p))])
upper_q = max(viz.df$log10p[!(is.na(viz.df$log10p))])
middle_q = mean(c(lower_q, upper_q)) # median(viz.df$log10p[!(is.na(viz.df$log10p))]) #
# middle_q = 7.5
g<-ggplot(data = viz.df, aes(x = Condition, y = Enrichment.Term, color = log10p, size = Frequency)) +
geom_point() +
# scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = '-log10(q-val)',
# limits = c(lower_q, upper_q), midpoint = middle_q) +
scale_colour_gradientn(colours = rev(c("darkred", "orange", "yellow")), name = '-log10(q-val)')+
scale_size_continuous(range = c(1,15)) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
text = element_text(size = 20))
g
Warning message: “Removed 141 rows containing missing values (`geom_point()`).”
First, we need to group genes into co-expression clusters. To do so, we will run hdWGCNA on the entire CD8+ dataset, then look at how the DE genes identified in the above analysis group by these co-expression clusters
Running on just the subset of DE genes or cell subsets in the branches is too underpowered for hdWGCNA
We run this on the raw counts and then normalize the metacells as suggested by the tutorial
wgcna.so <- abc.cd8s.all
DefaultAssay(wgcna.so)<-'RNA'
wgcna.so <- hdWGCNA::SetupForWGCNA(seurat_obj = wgcna.so,
gene_select = "fraction", # the gene selection approach
fraction = 0.05,
wgcna_name = 'CD8.wgcna')
# construct metacells in each group
wgcna.so <- MetacellsByGroups(
seurat_obj = wgcna.so,
group.by = c("Cell.Type.Level2", "orig.ident"), # specify the columns in seurat_obj@meta.data to group by
reduction = 'pca', # select the dimensionality reduction to perform KNN on
k = 20, # nearest-neighbors parameter
max_shared = 5, # maximum number of shared cells between two metacells
ident.group = 'Cell.Type.Level2', # set the Idents of the metacell seurat object
assay = 'RNA',
slot = 'counts'
)
wgcna.so <- NormalizeMetacells(wgcna.so)
# # normalize/batch correct the metacells
# suppressMessages({
# suppressWarnings({
# wgcna.so <- NormalizeMetacells(wgcna.so)
# wgcna.so <- ScaleMetacells(wgcna.so, features=VariableFeatures(abc.cd8s.all))
# wgcna.so <- RunPCAMetacells(wgcna.so, features=VariableFeatures(abc.cd8s.all))
# wgcna.so <- RunHarmonyMetacells(wgcna.so, group.by.vars='orig.ident')
# wgcna.so <- RunUMAPMetacells(wgcna.so, reduction='harmony', dims=1:15)
# })
# })
# p1 <- DimPlotMetacells(wgcna.so, group.by='Cell.Type.Level2') + umap_theme() + ggtitle("Cell Type")
# p2 <- DimPlotMetacells(wgcna.so, group.by='orig.ident') + umap_theme() + ggtitle("Sample")
# p1 | p2
wgcna.so <- hdWGCNA::SetDatExpr(seurat_obj = wgcna.so,
group_name = as.character(unique(wgcna.so$Cell.Type.Level2)),
group.by='Cell.Type.Level2',
use_metacells=TRUE,
assay = 'RNA', slot = 'data' # use the log-normalized counts
)
Warning message in MetacellsByGroups(seurat_obj = wgcna.so, group.by = c("Cell.Type.Level2", :
“Removing the following groups that did not meet min_cells: CD8+ TE/Ex#ABC, CD8+ TE/Ex#DT_Veh, CD8+ TE/Ex#UNTR, CD8+ TEA_1#DT_Veh, CD8+ TEA_1#UNTR, CD8+ TEA_2#ABC, CD8+ TEA_2#DT_ABC, CD8+ TEA_2#DT_Veh, CD8+ TEA_2#UNTR, CD8+ TEx prec#ABC, CD8+ TEx prec#DT_Veh, CD8+ TEx prec#UNTR, CD8+ TEx#ABC, CD8+ TEx#DT_Veh, CD8+ TEx#UNTR, CD8+ TN/EA-ISG#DT_Veh, CD8+ TN/EA-ISG#UNTR”
# Test different soft powers:
wgcna.so <- TestSoftPowers(
wgcna.so,
networkType = 'signed' # you can also use "unsigned" or "signed hybrid"
)
# plot the results:
plot_list <- hdWGCNA::PlotSoftPowers(wgcna.so)
# assemble with patchwork
patchwork::wrap_plots(plot_list, ncol=2)
pickSoftThreshold: will use block size 5219. pickSoftThreshold: calculating connectivity for given powers... ..working on genes 1 through 5219 of 8572 ..working on genes 5220 through 8572 of 8572 Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k. 1 1 0.34500 19.40 0.796 4.34e+03 4.34e+03 4620.00 2 2 0.00781 -1.44 0.944 2.22e+03 2.21e+03 2650.00 3 3 0.27200 -6.09 0.780 1.15e+03 1.13e+03 1600.00 4 4 0.71900 -7.22 0.862 6.02e+02 5.79e+02 1010.00 5 5 0.92600 -6.98 0.966 3.20e+02 3.01e+02 666.00 6 6 0.94800 -5.85 0.970 1.72e+02 1.57e+02 454.00 7 7 0.95600 -4.82 0.978 9.38e+01 8.28e+01 320.00 8 8 0.95700 -4.12 0.979 5.21e+01 4.38e+01 233.00 9 9 0.95800 -3.60 0.981 2.95e+01 2.33e+01 173.00 10 10 0.95100 -3.22 0.979 1.70e+01 1.25e+01 132.00 11 12 0.91100 -2.76 0.965 6.13e+00 3.65e+00 81.10 12 14 0.90700 -2.34 0.972 2.47e+00 1.10e+00 53.00 13 16 0.90000 -2.05 0.973 1.13e+00 3.37e-01 36.30 14 18 0.85700 -1.88 0.925 5.80e-01 1.05e-01 25.80 15 20 0.85700 -1.71 0.926 3.31e-01 3.37e-02 18.80 16 22 0.92500 -1.52 0.981 2.06e-01 1.10e-02 14.10 17 24 0.95800 -1.40 0.990 1.36e-01 3.67e-03 10.70 18 26 0.97300 -1.31 0.983 9.49e-02 1.25e-03 8.31 19 28 0.98000 -1.28 0.987 6.85e-02 4.30e-04 7.02 20 30 0.98700 -1.25 0.994 5.10e-02 1.51e-04 6.02 Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k. 1 1 0.344777420 19.393750 0.7955631 4336.2207 4337.3481 4624.2984 2 2 0.007808727 -1.436221 0.9442834 2219.6247 2205.2754 2653.7066 3 3 0.272279141 -6.086219 0.7800630 1149.3675 1126.7147 1602.5404 4 4 0.718665099 -7.217655 0.8624287 602.2537 579.3024 1012.2800 5 5 0.926057500 -6.976546 0.9658155 319.5523 300.7122 665.8585 6 6 0.947885285 -5.848773 0.9704821 171.8724 157.2604 454.2945
# construct co-expression network:
wgcna.so <- ConstructNetwork(
wgcna.so, soft_power=5,
setDatExpr=FALSE,
overwrite_tom = TRUE
)
modules <- GetModules(wgcna.so)
Calculating consensus modules and module eigengenes block-wise from all genes
Calculating topological overlaps block-wise from all genes
Flagging genes and samples with too many missing values...
..step 1
TOM calculation: adjacency..
..will use 20 parallel threads.
Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.
..Working on block 1 .
..Working on block 1 .
..merging consensus modules that are too close..
Number of genes assigned to each co-expression cluster (grey means they were not assigned)
table(modules$module)
grey green turquoise yellow brown blue red black
6160 172 703 197 530 547 110 94
pink
59
The distribution of these within the DE genes identified during the trajectory analysis:
gene.map<-modules$module
names(gene.map)<-rownames(modules)
de.res[['WGCNA.cluster']]<-gene.map[de.res$gene.id]
de.res[is.na(de.res$WGCNA.cluster), 'WGCNA.cluster'] <- 'grey'
table(de.res$WGCNA.cluster)
grey green turquoise yellow brown blue red black
503 0 11 3 18 33 27 6
pink
5
Compare this to all the DE genes:
table(de.res[c('Condition', 'WGCNA.cluster')])
WGCNA.cluster Condition grey green turquoise yellow brown blue red black pink Branch1.cells 78 0 1 1 3 4 5 1 1 Branch3.cells 137 0 3 0 3 7 6 1 1 Path1A.cells 128 0 1 0 3 5 5 2 1 Path1B.cells 118 0 4 1 3 7 6 1 1 Path2.cells 42 0 2 1 6 10 5 1 1
Next, let's see whether any of the conditions are enriched for any of the coexpression modules (disregarding "grey" genes not assigned to a pathway):
# p.adjust = BH corrected -- # https://github.com/YuLab-SMU/clusterProfiler/issues/124
# GeneRatio -- ratio of input genes in specific pathway vs input genes in all pathways https://www.biostars.org/p/220465/
# can improve visualization by directly manipulating the @result of each output
# if doing this, change qvalueCutoff to 1 and directly filter
visualize_ora<-function(cp.out, factor.name, sig.thresh = 0.1, top_n_terms = NULL){
viz.df<-cp.out@result
if (is.null(top_n_terms)){
top_n_terms = dim(viz.df)[[1]]
}
viz.df[['significant']]<-paste0('fdr > ', sig.thresh)
viz.df[viz.df$p.adjust <= sig.thresh, 'significant']<-paste0('fdr', ' \u2264 ', sig.thresh)
viz.df[['significant']]<-factor(viz.df[['significant']],
levels = c(paste0('fdr', ' \u2264 ', sig.thresh), paste0('fdr > ', sig.thresh)))
viz.df[['log.fdr']]<- -log10(viz.df$p.adjust)
viz.df[['GeneRatio2']]<-unlist(unname(sapply(viz.df$GeneRatio, function(x) eval(parse(text=x)))))
viz.df<-viz.df[with(viz.df, order(significant, -GeneRatio2, -log.fdr, -Count)), ]
viz.df<-viz.df[1:top_n_terms, ]
viz.df[['ID']]<-factor(viz.df[['ID']],
levels = rev(viz.df[['ID']]))
green_hex = brewer.pal(n = 11, name ='RdYlGn')[[11]]
yellow_hex = brewer.pal(n = 11, name = 'RdYlGn')[[6]]
red_hex = brewer.pal(n = 11, name = 'RdYlGn')[[1]]
lower_q = min(viz.df$log.fdr[!(is.na(viz.df$log.fdr))])
upper_q = max(viz.df$log.fdr[!(is.na(viz.df$log.fdr))])
middle_q = mean(c(lower_q, upper_q))#median(viz.df$log.fdr[!(is.na(viz.df$log.fdr))])
# middle_q = 7.5
g<-ggplot(data = viz.df, aes(x = GeneRatio2, y = ID, color = log.fdr, size = Count, shape = significant)) +
geom_point() +
ylab('') + xlab('GeneRatio') + ggtitle(factor.name) +
# scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = '-log10(fdr)',
# limits = c(lower_q, upper_q), midpoint = middle_q) +
scale_colour_gradientn(colours = rev(c("darkred", "orange", "yellow")), name = '-log10(fdr)')+
scale_shape_manual(values = c(16,4), drop = FALSE)+
scale_size_continuous(range = c(4,12))+
theme_bw() + theme(#axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
text = element_text(size = 15), plot.title = element_text(hjust = 0.5))
return(g)
}
visualize_ora_all<-function(ora_res, sig.thresh = 0.1, top_n_terms = 15) {
viz.df <- data.frame(matrix(ncol = 10))
colnames(viz.df)<-c('ID','Description','GeneRatio','BgRatio','pvalue','p.adjust','qvalue','geneID','Count', 'Comparison')
for (comp.name in names(ora_res))({
viz.df.sub<-ora_res[[comp.name]]@result
viz.df.sub[['Comparison']]<-comp.name
viz.df<-rbind(viz.df, viz.df.sub)
})
viz.df<-viz.df[-1,]
rownames(viz.df)<-1:dim(viz.df)[[1]]
viz.df[['significant']]<-paste0('fdr > ', sig.thresh)
viz.df[viz.df$p.adjust <= sig.thresh, 'significant']<-paste0('fdr', ' \u2264 ', sig.thresh)
viz.df[['significant']]<-factor(viz.df[['significant']],
levels = c(paste0('fdr', ' \u2264 ', sig.thresh), paste0('fdr > ', sig.thresh)))
viz.df[['log.fdr']]<- -log10(viz.df$p.adjust)
viz.df[['GeneRatio2']]<-unlist(unname(sapply(viz.df$GeneRatio, function(x) eval(parse(text=x)))))
viz.df<-viz.df[with(viz.df, order(significant, -GeneRatio2, -log.fdr, -Count)), ]
n_terms<-top_n_terms
unique_n_terms<-length(unique(viz.df[1:n_terms, 'ID']))
while ((unique_n_terms < top_n_terms) & (n_terms < dim(viz.df)[[1]])){
n_terms <- n_terms + 1
unique_n_terms<-length(unique(viz.df[1:n_terms, 'ID']))
}
viz.df<-viz.df[1:n_terms, ]
green_hex = brewer.pal(n = 11, name ='RdYlGn')[[11]]
yellow_hex = brewer.pal(n = 11, name = 'RdYlGn')[[6]]
red_hex = brewer.pal(n = 11, name = 'RdYlGn')[[1]]
lower_q = min(viz.df$log.fdr[!(is.na(viz.df$log.fdr))])
upper_q = max(viz.df$log.fdr[!(is.na(viz.df$log.fdr))])
middle_q = mean(c(lower_q, upper_q))#median(viz.df$log.fdr[!(is.na(viz.df$log.fdr))])
# middle_q = 7.5
g<-ggplot(data = viz.df, aes(x = Comparison, y = ID, color = log.fdr, size = GeneRatio2, shape = significant)) +
geom_point() +
ylab('') + xlab('') + ggtitle('') +
# scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = '-log10(fdr)',
# limits = c(lower_q, upper_q), midpoint = middle_q) +
scale_colour_gradientn(colours = rev(c("darkred", "orange", "yellow")), name = '-log10(fdr)')+
scale_shape_manual(values = c(16,4), drop = FALSE)+
scale_size_continuous(range = c(4,12))+
theme_bw() + theme(#axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
text = element_text(size = 20), plot.title = element_text(hjust = 0.5))
return(g)
}
term2gene <- modules[c('color', 'gene_name')]
colnames(term2gene) <- c('Pathway', 'gene.id')
term2gene <- term2gene[term2gene$Pathway != 'grey', ] # ignore the genes not assigned to a pathway
term2gene$Pathway <- factor(term2gene$Pathway)
rownames(term2gene) <- 1:dim(term2gene)[[1]]
We lower the minimum gene set size to 2 since there are so few genes:
ora_res <- list()
for (cond.name in unique(de.res$Condition)){
dr <- de.res[de.res$Condition == cond.name, ]
dr <- dr[dr$WGCNA.cluster != 'grey', ] # ignore genes not assigned to a pathway
dr <- dr$gene.id
ora_res[[cond.name]]<-clusterProfiler::enricher(dr, TERM2GENE = term2gene,
pAdjustMethod = 'BH',
pvalueCutoff = 1, # will filter later
qvalueCutoff = 1, # will filter later
minGSSize = 2,
maxGSSize = 1000
)
}
# visualize
dotplot_list<-c()
for (fn in names(ora_res)){
cp.out<-ora_res[[fn]]
g<-visualize_ora(cp.out, fn, top_n_terms = length(unique(term2gene$Pathway)))
dotplot_list[[fn]]<-g
}
h_ = 10
w_ = 35
options(repr.plot.height=h_, repr.plot.width=w_)
suppressWarnings({
g<-cowplot::plot_grid(dotplot_list[[1]], dotplot_list[[2]], dotplot_list[[3]], dotplot_list[[4]], dotplot_list[[5]],
ncol = 3)
})
g
h_ = 10
w_ = 20
options(repr.plot.height=h_, repr.plot.width=w_)
g<-visualize_ora_all(ora_res, sig.thres = 0.1, top_n_terms = 20)
g
Oddly enough, they are all enriched for the same co-expression cluster - the 'red' one
Next, let's visualize these co-expression clusters (median of the integraed, scaled counts across genes per cell) on the umap in the subregions for each trajectory analysis:
graph_branch_module<-function(cluster.label, ind.scale = T, scale = NULL,
de.genes, so, barcodes.branch, viz.df.0, assay, slot){
de.genes.cluster<-de.genes[de.genes$WGCNA.cluster == cluster.label, 'gene.id']
if (length(de.genes.cluster) > 0){
expr<-GetAssayData(object = so, slot = slot, assay = assay)[de.genes.cluster, barcodes.branch]
n.genes<-length(de.genes.cluster)
viz.df<-viz.df.0[barcodes.branch, ]
if (length(de.genes.cluster) > 1){
viz.df[['median.expression']]<-colMedians(expr)
}else{
viz.df[['median.expression']]<-unlist(unname(expr))
}
if (ind.scale){
lower_q = min(viz.df$median.expression[!(is.na(viz.df$median.expression))])
upper_q = max(viz.df$median.expression[!(is.na(viz.df$median.expression))])
middle_q = mean(c(lower_q, upper_q))
}else{
lower_q = scale[1]
upper_q = scale[2]
middle_q = scale[3]
}
g<-ggplot(viz.df, aes(x = UMAP_1, y = UMAP_2, color = median.expression)) + geom_point() +
scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = 'median expression',
limits = c(lower_q, upper_q), midpoint = middle_q)+
ggtitle(paste0('WGCNA Cluster: ', cluster.label, ' | no. genes: ', n.genes))+
theme_bw()
return(g)
}
else{
return(NULL)
}
}
slot = 'scale.data'
assay = 'integrated'
umap.coords<-abc.cd8s.all[["umap"]]@cell.embeddings
all.paths.viz<-list()
suppressWarnings({
for (branch.label in names(branch.cells.list)){
viz.df.0<-as.data.frame(umap.coords)
barcodes.branch<-branch.cells.list[[branch.label]]
de.genes<-de.res[de.res$Condition == branch.label, ]
# all genes
expr<-GetAssayData(object = abc.cd8s.all, slot = slot, assay = assay)[de.genes$gene.id, barcodes.branch]
n.genes<-dim(expr)[[1]]
viz.df<-viz.df.0
viz.df[['median.expression']]<-NA
viz.df[barcodes.branch, 'median.expression']<-colMedians(expr)
lower_q = min(viz.df$median.expression[!(is.na(viz.df$median.expression))])
upper_q = max(viz.df$median.expression[!(is.na(viz.df$median.expression))])
middle_q = mean(c(lower_q, upper_q))
g.all<-ggplot(viz.df, aes(x = UMAP_1, y = UMAP_2, color = median.expression)) + geom_point() +
scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = 'median expression',
limits = c(lower_q, upper_q), midpoint = middle_q)+
ggtitle(paste0(branch.label, ' - All DE genes', ' | no. genes: ', n.genes))+
theme_bw()
g.red<-graph_branch_module(cluster.label = 'red',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.blue<-graph_branch_module(cluster.label = 'blue',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.brown<-graph_branch_module(cluster.label = 'brown',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.pink<-graph_branch_module(cluster.label = 'pink',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.turquoise<-graph_branch_module(cluster.label = 'turquoise',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.black<-graph_branch_module(cluster.label = 'black',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.yellow<-graph_branch_module(cluster.label = 'yellow',
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
all.paths.viz[[branch.label]]<-cowplot::plot_grid(g.all, g.red, g.blue, g.brown, g.pink, g.turquoise, g.black, g.yellow,
ncol = 1, title = branch.label)
}
})
Each on different color scales:
h_ = 20
w_ = 30
options(repr.plot.height=h_, repr.plot.width=w_)
g.final<-cowplot::plot_grid(all.paths.viz[[1]], all.paths.viz[[2]], all.paths.viz[[3]], all.paths.viz[[4]],
all.paths.viz[[5]],
ncol = length(all.paths.viz))
g.final
slot = 'scale.data'
assay = 'integrated'
umap.coords<-abc.cd8s.all[["umap"]]@cell.embeddings
all.paths.viz<-list()
suppressWarnings({
for (branch.label in names(branch.cells.list)){
viz.df.0<-as.data.frame(umap.coords)
barcodes.branch<-branch.cells.list[[branch.label]]
de.genes<-de.res[de.res$Condition == branch.label, ]
# all genes
expr<-GetAssayData(object = abc.cd8s.all, slot = slot, assay = assay)[de.genes$gene.id, barcodes.branch]
n.genes<-dim(expr)[[1]]
viz.df<-viz.df.0
viz.df[['median.expression']]<-NA
viz.df[barcodes.branch, 'median.expression']<-colMedians(expr)
lower_q = -2
upper_q = 2
middle_q = 0
g.all<-ggplot(viz.df, aes(x = UMAP_1, y = UMAP_2, color = median.expression)) + geom_point() +
scale_color_gradient2(low = green_hex, mid = yellow_hex, high = red_hex, name = 'median expression',
limits = c(lower_q, upper_q), midpoint = middle_q)+
ggtitle(paste0(branch.label, ' - All DE genes', ' | no. genes: ', n.genes))+
theme_bw()
g.red<-graph_branch_module(cluster.label = 'red', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.blue<-graph_branch_module(cluster.label = 'blue', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.brown<-graph_branch_module(cluster.label = 'brown', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.pink<-graph_branch_module(cluster.label = 'pink', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.turquoise<-graph_branch_module(cluster.label = 'turquoise', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.black<-graph_branch_module(cluster.label = 'black', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
g.yellow<-graph_branch_module(cluster.label = 'yellow', ind.scale = F, scale = c(lower_q, upper_q, middle_q),
de.genes = de.genes, so = abc.cd8s.all, barcodes.branch = barcodes.branch, viz.df.0 = viz.df.0, assay = assay, slot = slot)
all.paths.viz[[branch.label]]<-cowplot::plot_grid(g.all, g.red, g.blue, g.brown, g.pink, g.turquoise, g.black, g.yellow,
ncol = 1, title = branch.label)
}
})
All on the same color scale:
h_ = 20
w_ = 30
options(repr.plot.height=h_, repr.plot.width=w_)
g.final<-cowplot::plot_grid(all.paths.viz[[1]], all.paths.viz[[2]], all.paths.viz[[3]], all.paths.viz[[4]],
all.paths.viz[[5]],
ncol = length(all.paths.viz))
g.final